%% ====================== [Monte Carlo Method for Distribution] ==========================
clear
vMat=[32.5e-3,pi*0.1^2,0.467,3];
coeOfVar=[0 0 0 0.2];
svm=4;
usvm=[1,2,3];
vNum=length(svm);
simMCNum=500;
disp('Pavement sine signal - input 1')
disp('Pavement random signal - input 2')
exc=input('signal:');
swiMat=[-1,1];
excSwitch=swiMat(exc);
%----------------------------- [model parameter area] -------------------------------
Pa=101*10^3;
A_level=16; %-|
B_level=64; %-|-pavement
C_level=256; %-|
Level=B_level;
v=60; %Vehicle speed
C=2000; %damping   
m=800;
randomNum=simMCNum^vNum;
accMCData=zeros(10001,randomNum);
devMCData=zeros(10001,randomNum);
%--------------------------------------------------------------------------
for index=1:randomNum
    [d0,Ac,L0,k]=MCiniRandom(vMat,coeOfVar,svm,...
    usvm,index,simMCNum);
    sim('quasi0StiffSuspension.slx')
    accMCData(:,index)=acc.data;
    devMCData(:,index)=dev.data;
    %----------------------------------------------------------------------
    proPers=floor(index/randomNum*100); %
    clc %---progress bar
    disp(['[','*'*ones(1,proPers-1),'>',' '*ones(1,100-proPers) ...       %
        ,']',num2str(proPers),'%']) %
    % ----------------------------------------------------------------------
end
%% ======================== [a certain second probability distribution] =========================
[Xtpdf,Xt] = ksdensity(devMCData(3501,:));
figure(6)
plot(Xt,Xtpdf,'LineWidth',1)
title('Output probability density at 5s')
xlabel('output random variable')
ylabel('probability density \it{f}')
%--------------------------------------------------------------------------
clc
disp(['[','<'*ones(1,proPers),' '*ones(1,100-proPers) ...             %
        ,']',num2str(proPers),'%'])